library(gplots)

params <- matrix(nrow=1080,ncol=6)
dpl <- c(0,0.25,0.5,0.75,1)
scl <- c(0.001,0.01,0.1,0.5,1,2)
sml <- c(10^-5,10^-6,10^-7,10^-8)
nml <- c(10^-4,10^-5,10^-7)
sw <- 0

for (divprob in 1:5){
  for (selcoeff in 1:6){
    for (selmu in 1:4){
      for (neutmu in 1:3){
        for (replicates in 1:3){
          sw <- sw + 1
          params[sw,1] <- sw
          params[sw,2] <- dpl[divprob]
          params[sw,3] <- scl[selcoeff]
          params[sw,4] <- sml[selmu]
          params[sw,5] <- nml[neutmu]
          params[sw,6] <- replicates
        }
      }
    }	
  }
}

smu <- 10^-7
# Get the three panels of s=0.5 s=1, s=2
s <- 0.5
res1 <- matrix(data=NA,nrow=8,ncol=45)
column <- 0
for (r in 1:length(dpl)){
  for (nmu in 1:length(nml)){
    for (repl in 1:3){
      column <- column+1
      sw <- params[params[,3]==s&params[,4]==smu&params[,2]==dpl[r]&params[,5]==nml[nmu]&params[,6]==repl,1]
      f <- read.csv(paste("../sim",sw,"/sim",sw,".hit.log",sep=""),header=FALSE)
      f <- f[f[,1]==7300,]
      res1[1:5,column] <- as.numeric(f[,3:7])
      # get the selective clones
      d <- read.csv(paste("../sim",sw,"/sim",sw,".frequencies",sep=""),header=FALSE)
      # get clone_ids and frequencies that satisfy criteria: end of simulation, frequency higher than 1%, selective clone (last hit is a selective locus)
      dat <- d[d[,2]==7300&d[,3]>0.01&d[,1]>0&d[,4]>104&d[,4]<110,c(1,3)]
      # get s from clone_ids, from .clones the number of beneficial hits that the clone has
      bh <- read.csv(paste("../sim",sw,"/sim",sw,".clones",sep=""),header=FALSE)
      # match subset to superset, get indices of the subset elements in the superset
      ind <- match(dat[,1],bh[,1])
      clones.s <- bh[ind,c(1,3)]

      # Get index of clones with maximum s -> last category, the most hits
      clones.s <- clones.s[which(clones.s[,2]==max(clones.s[,2],na.rm=TRUE)),]
      #if(nrow(clones.s)<1){print (paste("this:",sw,r,nmu,repl))}
      number.interfering.clones <- nrow(clones.s)
      ind <- match(clones.s[,1],dat[,1])
      avg.freq <- mean(dat[ind,2])
      sd.freq <- sd(dat[ind,2])
      res1[6,column] <- number.interfering.clones
      res1[7,column] <- avg.freq
      res1[8,column] <- sd.freq     
    }
  }
}

# Get the three panels of s=0.5 s=1, s=2
s <- 1
res2 <- matrix(data=NA,nrow=8,ncol=45)
column <- 0
for (r in 1:length(dpl)){
  for (nmu in 1:length(nml)){
    for (repl in 1:3){
      column <- column+1
      sw <- params[params[,3]==s&params[,4]==smu&params[,2]==dpl[r]&params[,5]==nml[nmu]&params[,6]==repl,1]
      f <- read.csv(paste("../sim",sw,"/sim",sw,".hit.log",sep=""),header=FALSE)
      f <- f[f[,1]==7300,]
      res2[1:5,column] <- as.numeric(f[,3:7])
      # get the selective clones
      d <- read.csv(paste("../sim",sw,"/sim",sw,".frequencies",sep=""),header=FALSE)
      # get clone_ids and frequencies that satisfy criteria: end of simulation, frequency higher than 1%, selective clone (last hit is a selective locus)
      dat <- d[d[,2]==7300&d[,3]>0.01&d[,1]>0&d[,4]>104&d[,4]<110,c(1,3)]
      # get s from clone_ids, from .clones the number of beneficial hits that the clone has
      bh <- read.csv(paste("../sim",sw,"/sim",sw,".clones",sep=""),header=FALSE)
      # match subset to superset, get indices of the subset elements in the superset
      ind <- match(dat[,1],bh[,1])
      clones.s <- bh[ind,c(1,3)]

      # Get index of clones with maximum s -> last category, the most hits
      clones.s <- clones.s[which(clones.s[,2]==max(clones.s[,2],na.rm=TRUE)),]
      #if(nrow(clones.s)<1){print (paste("this:",sw,r,nmu,repl))}
      number.interfering.clones <- nrow(clones.s)
      ind <- match(clones.s[,1],dat[,1])
      avg.freq <- mean(dat[ind,2])
      sd.freq <- sd(dat[ind,2])
      res2[6,column] <- number.interfering.clones
      res2[7,column] <- avg.freq
      res2[8,column] <- sd.freq     
    }
  }
}

s <- 2
res3 <- matrix(data=NA,nrow=8,ncol=45)
column <- 0
for (r in 1:length(dpl)){
  for (nmu in 1:length(nml)){
    for (repl in 1:3){
      column <- column+1
      sw <- params[params[,3]==s&params[,4]==smu&params[,2]==dpl[r]&params[,5]==nml[nmu]&params[,6]==repl,1]
      f <- read.csv(paste("../sim",sw,"/sim",sw,".hit.log",sep=""),header=FALSE)
      f <- f[f[,1]==7300,]
      res3[1:5,column] <- as.numeric(f[,3:7])
      # get the selective clones
      d <- read.csv(paste("../sim",sw,"/sim",sw,".frequencies",sep=""),header=FALSE)
      # get clone_ids and frequencies that satisfy criteria: end of simulation, frequency higher than 1%, selective clone (last hit is a selective locus)
      dat <- d[d[,2]==7300&d[,3]>0.01&d[,1]>0&d[,4]>104&d[,4]<110,c(1,3)]
      # get s from clone_ids, from .clones the number of beneficial hits that the clone has
      bh <- read.csv(paste("../sim",sw,"/sim",sw,".clones",sep=""),header=FALSE)
      # match subset to superset, get indices of the subset elements in the superset
      ind <- match(dat[,1],bh[,1])
      clones.s <- bh[ind,c(1,3)]

      # Get index of clones with maximum s -> last category, the most hits
      clones.s <- clones.s[which(clones.s[,2]==max(clones.s[,2],na.rm=TRUE)),]
      #if(nrow(clones.s)<1){print (paste("this:",sw,r,nmu,repl))}
      number.interfering.clones <- nrow(clones.s)
      ind <- match(clones.s[,1],dat[,1])
      avg.freq <- mean(dat[ind,2])
      sd.freq <- sd(dat[ind,2])
      res3[6,column] <- number.interfering.clones
      res3[7,column] <- avg.freq
      res3[8,column] <- sd.freq     
    }
  }
}

write.table(res1,"res1.sm10-5.csv",sep=",",quote=FALSE,row.names=FALSE,col.names=FALSE)
write.table(res2,"res2.sm10-5.csv",sep=",",quote=FALSE,row.names=FALSE,col.names=FALSE)
write.table(res3,"res3.sm10-5.csv",sep=",",quote=FALSE,row.names=FALSE,col.names=FALSE)

write.table(res1,"res1.sm10-6.csv",sep=",",quote=FALSE,row.names=FALSE,col.names=FALSE)
write.table(res2,"res2.sm10-6.csv",sep=",",quote=FALSE,row.names=FALSE,col.names=FALSE)
write.table(res3,"res3.sm10-6.csv",sep=",",quote=FALSE,row.names=FALSE,col.names=FALSE)

write.table(res1,"res1.sm10-7.csv",sep=",",quote=FALSE,row.names=FALSE,col.names=FALSE)
write.table(res2,"res2.sm10-7.csv",sep=",",quote=FALSE,row.names=FALSE,col.names=FALSE)
write.table(res3,"res3.sm10-7.csv",sep=",",quote=FALSE,row.names=FALSE,col.names=FALSE)


pdf("beneficial.mutations.frequencies.sm10-7.pdf",width=6.5, height=3)
par(mfrow=c(4,1))
for(i in 1:3){
  if(i==1){
    res <- res1
    par(mar=c(0,4,0.5,1))}
  if(i==2){
    res <- res2
    par(mar=c(0.5,4,0.5,1))}
  if(i==3){
    res <- res3
    par(mar=c(0,4,0.5,1))}
subres <- res[,1:9]
subres1 <- subres[,order(subres[1,]*1+subres[2,]*2+subres[3,]*4+subres[4,]*8+subres[5,]*16,decreasing=FALSE)]
subres <- res[,10:18]
subres2 <- subres[,order(subres[1,]*1+subres[2,]*2+subres[3,]*4+subres[4,]*8+subres[5,]*16,decreasing=FALSE)]
subres <- res[,19:27]
subres3 <- subres[,order(subres[1,]*1+subres[2,]*2+subres[3,]*4+subres[4,]*8+subres[5,]*16,decreasing=FALSE)]
subres <- res[,28:36]
subres4 <- subres[,order(subres[1,]*1+subres[2,]*2+subres[3,]*4+subres[4,]*8+subres[5,]*16,decreasing=FALSE)]
subres <- res[,37:45]
subres5 <- subres[,order(subres[1,]*1+subres[2,]*2+subres[3,]*4+subres[4,]*8+subres[5,]*16,decreasing=FALSE)]
nc <- t(t(c(NA,NA,NA,NA,NA,NA,NA,NA)))
ress <- cbind(subres1,nc,subres2,nc,subres3,nc,subres4,nc,subres5)
scheme <- c("magenta","blue","green","orange","red")
ress <- as.matrix(ress)
# split matrix, remove number of interfering clones
ci.labels <- ress[6:8,]
ress <- ress[1:5,]
rownames(ress) <- c("1 hit","2 hits","3 hits","4 hits", "5 hits")
mp <- barplot2(ress, beside = FALSE,col = scheme[1:5],ylim = c(0,1), plot.grid = FALSE,axes=FALSE)
axis(2,at=c(0,.5,1),label=c(0,.5,1))
#if(i==3){axis(1,at=c(c(5,15,25,35,45)*1.2-0.5),label=c("r=0","r=.25","r=.5","r=.75","r=1"),tick=FALSE)}
# calculate y-values
y.value.labels <- apply(ress[1:5,],2,sum)
#text(x=c(1:49)*1.2-0.5,y=y.value.labels+0.05,ci.labels[1,1:45],cex=0.6)
#text(x=c(1:45)*1.2-0.5,y=y.value.labels+0.2,format(ci.labels[2,1:45],digits=1),cex=0.6)
#text(x=c(1:45)*1.2-0.5,y=y.value.labels+0.1,format(ci.labels[3,1:45],digits=1),cex=0.6)
#legend(24,1,legend=rownames(res),fill=gray(1-(0:4/4)))
  if(i==1){mtext(side=2, text="s=0.5",line=2)}
  if(i==2){mtext(side=2, text="s=1",line=2)}
  if(i==3){mtext(side=2, text="s=2",line=2)}
}
par(mar=c(0,4,0,1))
mp <- barplot2(ress, beside = FALSE,col = scheme[1:5],ylim = c(0,1), plot.grid = FALSE,axes=FALSE,plot=FALSE)
axis(1,at=c(c(5,15,25,35,45)*1.2-0.5),label=c("r=0","r=.25","r=.5","r=.75","r=1"),tick=FALSE)
#mtext(side=1, text="Variation between individual runs",line=2)
dev.off()


